Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des opérations statistiques avancées grâce au langage R. Ainsi, notre choix s’est porté sur des données statistiques relatives au prestige (score attribuée) de différentes professions canadiennes.
Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer une modélisation statistique sur notre jeu de données, d’etudier la coréalation entre des variables en effectuant une ou plusieurs regréssion linéaire et produire un rapport.
R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.
La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.
L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.
Notre jeu de données comporte 102 individus décrits par 6 variables :
head(data) %>% paged_table()viz_edu <- data %>%
ggplot(aes(x = prestige, y = education, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(6, 16, 2)
)
viz_edu <- ggplotly(viz_edu)
viz_inc <- data %>%
ggplot(aes(x = prestige, y = income, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 20000, 10000)
)
viz_inc <- ggplotly(viz_inc)
viz_women <- data %>%
ggplot(aes(x = prestige, y = women, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 100, 25)
)
viz_women <- ggplotly(viz_women)
viz_census <- data %>%
ggplot(aes(x = prestige, y = census, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(2500, 7500, 2500)
)
viz_census <- ggplotly(viz_census)
subplot(viz_edu, viz_inc, viz_women, viz_census, nrows = 2, margin = 0.07)Dans un premier temps, nous pouvons observer une forte relation linéaire entre prestige et les variables education et income, contrairement à women et census.
Dans un second temps, nous pouvons étudier d’éventuelles corrélations entre deux variables à l’aide d’un nuage de points. Ici, nous utilisons la fonction ggpairs afin de pouvoir en observer plusieurs en même temps.
generate_scatterM <- function(data, mapping){
data %>% ggplot(mapping = mapping) + geom_point(size = 0.8) +
geom_smooth(method = loess, color = "red", size = 0.85, se = FALSE) +
geom_smooth(method = lm, color = "blue", size = 0.85, se = FALSE)
}
scatterM <- data %>% ggpairs(columns = 1:5, lower = list(continuous = generate_scatterM))
scatterM <- ggplotly(scatterM)
scatterMeducation <- data %>%
ggplot(aes(x = education, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(20, 80, 20)
)
scatter_edu <- ggplotly(education)
scatter_eduLa droite bleue est définie par la méthode des moindres carrés (MSE), il s’agit d’une droite de régression linéaire entre les variables education et prestige.
La courbe rouge indique la tendance globale entre ces deux variables, il s’agit d’une courbe de régression de type lowess. Les deux courbes extérieures représentent l’intervalle de confiance de cette courbe de régression.
On constate que que la droite de régression est presque toujours comprise dans l’intervalle de confiance, l’hypothèse de linéarite entre les variables education et prestige est donc acceptable.
Nous testons aussi la relation entre les variables income et prestige.
income <- data %>%
ggplot(aes(x = income, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(0, 25000, 5000)
) +
scale_y_continuous(
breaks = seq(20, 100, 20)
)
scatter_inc <- ggplotly(income)
scatter_incIci, en regardant la forme du lien entre la variable entre les deux variables, on s’aperçoit que la droite de régression suis bien moins l’intervalle de confiance de la courbe lowess. l’hypothèse de linéarité est alors plus critiquable.
Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction “lm”.
lin_model_edu <- lm(prestige ~ education, data = data)
lin_model_edu##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Coefficients:
## (Intercept) education
## -10.732 5.361« Intercept » correspond ici à l’ordonnée à l’origine, le « b » de notre droite, et le « x » est la pente de la droite, ce qui correspond au « b » dans notre notation.
L’équation, de notre droite est donc \(y = -10.732 + 5.361x\).
Il existe différentes manières, entre diagrammes et tests statistiques, visant à évaluer le lien linéaire entre deux variables. Ce lien est dit « significatif » s’il remplit certaines conditions. En effet, les résidus doivent être indépendants, distribués selon une loin Normale de moyenne 0 et de façon homogène (variance constante).
En premier lieu, il faut évaluer l’hypothèse d’indépendance des résidus. À l’aide d’un lag plot, il possible de mettre en évidence la présence d’une auto-corrélation.
bacf <- acf(residuals(lin_model_edu), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(lin_model_edu)))
lag_plot <- bacfdf %>%
ggplot(aes(x = lag, y = acf)) +
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0)) +
geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +
geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
theme_bw()
lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plotSur le graphique ci-dessus, chaque trait correspond à un lag, ou coefficient de corrélation entre les résidus de chaque point, les pointillés bleus quant à eux sont les intervalles de confiance du coefficient de corrélation égale à 0.
En s’appuyant sur ce graphique, nous pouvons observer qu’une auto-corrélation significative est présente pour les lags 1 à 3, 5, 7 et 12 par exemple.
Afin d’aller plus loin et de détecter une éventuelle auto-corrélation des erreurs d’ordre 1 (lag = 1), il est possible d’emloyer le test de Durbin-Watson.
L’intérêt d’avoir recours à ce test est de déterminer si le modèle est perfectible. Souvent, l’autocorrélation s’observe sur les résidus d’une modélisation de série chronologique. Toutefois, si le type de modèle est bien choisi (en l’occurrence une régrassion linéaire simple), il peut exister une véritable autocorrélation entre les observations. En revanche, si aucune autocorrélation n’a de raison d’être, il faut chercher d’autres variables candidates.
durbinWatsonTest(lin_model_edu)## lag Autocorrelation D-W Statistic p-value
## 1 0.2752512 1.43917 0
## Alternative hypothesis: rho != 0Le test de Durbin-Watson indique qu’il existe une auto-corrélation significative entre les résidus d’une ligne du tableau de données et ceux de la ligne suivante.
La deuxième condition susmentionnée concerne la distribution des résidus selon une loi Normale de moyenne 0. La pertinence de l’ajustement d’une distribution donnée à un modèle peut être représentée graphiquement grâce à un diagramme quantile-quantile, ou QQ plot.
qq_edu <- lin_model_edu %>%
ggplot(aes(sample = rstandard(.))) +
stat_qq_line(color = "red", linetype = "dashed") +
stat_qq(size = 1.2) +
theme_bw() +
scale_x_continuous(
breaks = seq(-2, 2, 1)
) +
scale_y_continuous(
breaks = seq(-3, 3, 1)
) +
labs(title = "Normal Q-Q",
x = "Theoretical Quantiles\nlm(prestige ~ education)",
y = "Standardized residuals")
qq_edu <- ggplotly(qq_edu)
qq_eduNous constatons que les résidus sont relativement bien alignés le long de la droite figurant sur le graphique, nous pouvons donc en conclure que la distribution de cette série de données suit une loi normale.
En statistiques, le test de Shapiro–Wilk teste l’hypothèse nulle selon laquelle un échantillon de données est issu d’une population normalement distribuée. Ce test peut être employé afin d’évaluer l’hypothèse de normalité des résidus.
shapiro.test(residuals(lin_model_edu))##
## Shapiro-Wilk normality test
##
## data: residuals(lin_model_edu)
## W = 0.98065, p-value = 0.1406On considère que cette hypothèse est rejetée si la p-valeur est inférieure à \(0,05\), mais dans le cas présent le test accepte la normalité.
Il est essentiel d’appliquer des tests statistiques tels que celui de Shapiro-Wilk. Un exemple parlant serait celui d’une régression linéaire simple entres les variables income et prestige.
lin_model_inc <- lm(prestige~income, data = data)
lin_model_inc##
## Call:
## lm(formula = prestige ~ income, data = data)
##
## Coefficients:
## (Intercept) income
## 27.141176 0.002897qq_inc <- lin_model_inc %>%
ggplot(aes(sample = rstandard(.)))+
stat_qq_line(color = "red", linetype = "dashed") +
stat_qq(size = 1.2) +
theme_bw() +
scale_x_continuous(
breaks = seq(-2, 2, 1)
) +
scale_y_continuous(
breaks = seq(-3, 3, 1)
) +
labs(title = "Normal Q-Q",
x = "Theoretical Quantiles\nlm(prestige ~ income)",
y = "Standardized residuals")
qq_inc <- ggplotly(qq_inc)
qq_incL’alignement des résidus est légèrement moins précis que sur le QQ plot entre les variables education et prestige, mais il semble plutôt satisfaisant. On ne constaste aucune valeur aberrante.
shapiro.test(residuals(lin_model_inc))##
## Shapiro-Wilk normality test
##
## data: residuals(lin_model_inc)
## W = 0.97169, p-value = 0.02729Cependant, le test de Shapiro-Wilk renvoie une p-valeur inférieure à \(0,05\) et rejette donc l’hypothèse de normalité.
Enfin, le concept d’homoscédasticité est utilisé dans le contexte de la régression linéaire pour décrire le cas où la variance des erreurs du modèle est la même pour toutes les observations. Autrement dit, les variances sont homogènes et les erreurs identiquement distribuées. Il s’agit de l’une des propriétés fondamentales du modèle de régression linéaire.
Afin de visualiser cela, nous avons besoin des valeurs prédites par le modèle (fitted values).
Dans notre cas, le graphique représentera les valeurs de prestige prédite par le modèle pour les valeurs de education présentes dans les données.
sc_loc_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), sqrt(abs(rstandard(.))))) +
geom_point(size = 1.2) +
geom_smooth(method = loess, se = FALSE) +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(0, 1.5),
breaks = seq(0, 1.5, 0.5)
) +
labs(title = "Scale-location",
x = "Fitted values\nlm(prestige ~ education)",
y = "sqrt |Standardized residuals|")
sc_loc_edu <- ggplotly(sc_loc_edu)
sc_loc_eduLe graphique nous montre que les résidus sont répartis de façon homogène le long du gradient des valeurs prédites. Ceci est mis en évidence par la courbe de régression locale qui est quasiment plate.
Le test de Breusch-Pagan permet de tester l’hypothèse d’homoscédasticité du terme d’erreur d’un modèle de régression linéaire. Si la variance est constante, alors on a de l’homoscédasticité ; en revanche, si elle va varie, on a de l’hétéroscédasticité.
ncvTest(lin_model_edu)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6327545, Df = 1, p = 0.42635La p-valeur associée au test est supérieure au seuil de \(0.05\), ainsi l’hypothèse d’homoscédasticité est acceptée.
Afin de vérifier l’hypothèse selon laquelle la relation entre les variables education et prestige est linéaire, ainsi que l’homoscédasticité, il est possible de générer un residuals vs. fitted plot.
res_v_fit_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), residuals(.))) +
geom_point() +
stat_smooth(method="loess", se = FALSE) +
geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(-30, 20),
breaks = seq(-30, 20, 10)
) +
labs(title = "Residuals vs. Fitted values",
x = "Fitted values\nlm(prestige ~ education)",
y = "Residuals")
res_v_fit_edu <- ggplotly(res_v_fit_edu)
res_v_fit_eduIl n’y a pas de schéma clair dans les données, ni d’aberrations évidentes. Les résidus sont, à peu de choses près, uniformément distribués le long de la ligne 0. La droite de régression est bien adaptée aux données, par conséquent l’hypothèse de linéarité est acceptable.
En s’appuyant sur les différentes représentations graphiques ainsi que les tests effectués précédémment, on peut en conclure que les propriétés fondamentales de la régression linéaire (indépendance des résidus, normalité, homoscédasticité et linéarité) sont satisfaites.
Partant de ce postulat, les résultats de la régression linéaire entre les variables education et prestige peuvent être interprétés.
summary(lin_model_edu)##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.0397 -6.5228 0.6611 6.7430 18.1636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.732 3.677 -2.919 0.00434 **
## education 5.361 0.332 16.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
## F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16La fonction Summary nous permet d’avoir un bon apercus de ses resultats. Nous avons dans la première partie nomée Residuals, la norme des residus (min, max, quartile et mediane). La seconde partie nomée Coefficients est un tableau à deux ligne. La première ligne concerne l’ordonnée à l’origine, alors que la seconde concerne la pente. Le plus important dans notre cas est le coefficient de la pente égal à 5.361. Cela signifie que lorsque le niveau d’éducation augmente d’une unité, alors, le niveau du score de prestige augmente de 5.361 unités.
La fonction “confint” permet d’obtenir facilement l’intervalles de confiance des coefficients des paramètres. Cette fonction est basés sur la loi de Student.
confint(lin_model_edu) %>% kable()| 2.5 % | 97.5 % | |
|---|---|---|
| (Intercept) | -18.027220 | -3.436744 |
| education | 4.702223 | 6.019533 |
| l’intervalle d | e confiance à | 97.5% de la pente est une étendue de valeurs qui a une probabilité de 97.5% de contenir la vraie pente. |
On peut donc maintenant obtenir la réponse prédite par le modèle de régression, pour une valeur spécifique de la variable prédictive, qui n’a pas été observée. Par exemple, imaginons que je souhaite obtenir le score de prestige prédit par le modèle, pour un niveau d’éducation de 10.25, qui n’a pas été observé. Pour cela, il est nécessaire de créer un data frame contenant une variable éducation, avec les valeurs souhaitées. Ce nouveau data frame est ensuite passé dans l’argument “newdata” de la fonction “predict”.
predictions <- tibble(education=c(9.21, 11.07, 14.62))
predict(lin_model_edu, newdata = predictions, interval = "confidence") %>% kable()| fit | lwr | upr |
|---|---|---|
| 38.64170 | 36.58966 | 40.69374 |
| 48.61293 | 46.81134 | 50.41452 |
| 67.64405 | 64.52387 | 70.76423 |
Comme déjà vu, précédemment, les résidus de la régression peuvent être obtenus avec la fonction “‘residuals”.
data$residuals <- residuals(lin_model_edu)
head(data) %>% kable()| education | income | women | prestige | census | type | residuals |
|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 |
Les fitted, correspondent, aux prédictions du modèle de régression, mais pour les valeurs observées de la variable prédictive :
data$fitted <- fitted(lin_model_edu)
head(data) %>% kable()| education | income | women | prestige | census | type | residuals | fitted |
|---|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 | 59.54913 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 | 54.99238 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 | 57.72643 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 | 50.48924 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 | 67.64405 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 | 73.11215 |
La matrice de variance-covariance des paramètres du modèle peut s’obtenir avec la fonction “vcov” :
vcov(lin_model_edu) %>% kable()| (Intercept) | education | |
|---|---|---|
| (Intercept) | 13.520978 | -1.1835055 |
| education | -1.183506 | 0.1102162 |
pred_interval <- predict(lin_model_edu, interval = "prediction")
data <- cbind(data, pred_interval)linear_regression <- data %>%
ggplot(aes(y = prestige, x = education)) +
geom_point(size = 1.2) +
geom_smooth(color = "red", method = "lm", fill = "red") +
geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
geom_line(aes(y = upr), color = "red", linetype = "dashed") +
annotate("text", x = 9, y = 80,
label = paste0("prestige = ", coef(lin_model_edu)["(Intercept)"] %>% round(3),
" + ", coef(lin_model_edu)["education"] %>% round(3), " * education")) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(25, 75, 25)
) +
labs(title = "Linear Regression",
x = "Education", y = "Prestige")
linear_regression <- ggplotly(linear_regression)
linear_regressionLe modèle de régression linéaire est aussi bien utilisé pour chercher à prédire un phénomène que pour chercher à l’expliquer. Après avoir estimé un modèle de régression linéaire, on peut prédire quel serait le niveau de y pour des valeurs particulières de x. Il permet également d’estimer l’effet d’une ou plusieurs variables sur une autre en contrôlant par un ensemble de facteurs. En apprentissage statistique, la méthode de régression linéaire est considérée comme une méthode d’apprentissage supervisé utilisée pour prédire une variable quantitative. Dans cette perspective, on entraîne généralement le modèle sur un échantillon d’apprentissage et on teste ensuite les performances prédictives du modèle sur un échantillon de test.
Ce projet a été réalisé dans le cadre du cours de statistiques appliquées. Celui-ci convoque des connaissances en statistiques acquise plus tot cette année et fait le lien avec notre module d’apprentissage supervisé. En effet afin de concevoir nos systèmes prévisionnels, comprendre comment les relations entre les variables d’un jeu de donnée peuvent être étudier est une étape primordial. Il a été l’occasion pour nous, étudiants, de consolider nos compétences en analyse et manipulation de données, puis d’en acquérir de nouvelles par le biais de ce projet.